Chapter 12 High-level functions
In most situations you will have rasters that are larger than what you need. You might be interested in the landcover of your study area but the data covers all of Europe. Conversely, you might also need to combine different rasters if the scale of your analysis is larger than the data you are using. Often this is the case when larger rasters are split in smaller areas to reduce their size.
In the first case you would use crop().
Provided with a raster and an extend, crop() will return only the part of the raster that lies within the extend.
An extend usually is a named vector with four elements: xmin, xmax, ymin and ymax.
It might for example be the bounding box you derived from the vector data you were analyzing in sf.
You can create this raster yourself if you know, or can derive, the relevant numbers.
An alternative is the drawExtent() function from the raster package.
This lets you click on an already plotted raster an returns the extend of the thus created rectangle.
Try the code below in your own R instance to see what I mean.
plot(dem1)
extent = raster::drawExtent()
dem1crop = terra::crop(x =dem1, y=extent)

Nice!
Now have reduced the raster to the area we are actually interested in.
We can compute the mean hight or the sum of all hights with global().
global(dem1crop, "sum")## sum
## DGM25_2905530 100724.5
global(dem1crop, "mean")## mean
## DGM25_2905530 331.3307
We can also expand a raster with the extend() function.
That means we add rows and/ or columns with NA values.
The other way around we can also trim them, removing all cells with NAs in the margins.
dem1.exp = extend(dem1, y = c(100,100))
plot(dem1.exp)
Since all the NAs are just white, it’s hard to see what we actually did here.
Let’s give them some random value so we can see the cells we added.
na_cells = which(is.na(values(dem1.exp)))
values(dem1.exp)[na_cells] = 300
plot(dem1.exp)
We can now clearly see the added cells in brown.
With trim() we can remove them again.
na_cells = which(values(dem1.exp) == 300)
values(dem1.exp)[na_cells] = NA
dem1.trm = trim(dem1.exp)
plot(dem1.trm)
As mentioned before, sometimes you want to combine rasters.
For example our two DEMs are both part of a larger DEM that covers all of the German federal state Rhineland-Palatinate.
We can combine the two rasters with merge()
dem_merge = merge(dem1, dem2)
plot(dem_merge)
Whenever you want to use multiple raster in an operation, it is important that they have the same CRS and cell sizes.
You can increase cell sizes with aggregate() and reduce cell sizes with disaggregate().
In the first case, several cells are combined into one.
The number of cells is set with the fact argument.
When you supply one number it is used for both directions.
So if you call aggregate(x, fact = 2) and raster has the cell size 10m by 10m, the new raster will be 20m by 20m.
You can also provide fact with a vector that specifies the aggregation factor in each direction.
The third argument is fun which specifies how to calculate the new cell value.
If our new cell consists of four old cells (i.e., we used an aggregation factor of 2) we somehow need to combine those four numbers into one.
Common options here are the mean, the median, the modal, the maximum, the minimum, or the sum.
dem1_agg_2_mean = aggregate(dem1, fact=2, fun = "mean")
dem1_agg_3_mean = aggregate(dem1, fact=4, fun = "mean")
dem1_agg_4_mean = aggregate(dem1, fact=6, fun = "mean")
dem1_agg_5_mean = aggregate(dem1, fact=8, fun = "mean")
par(mfrow = c(2,2))
plot(dem1_agg_2_mean)
plot(dem1_agg_3_mean)
plot(dem1_agg_4_mean)
plot(dem1_agg_5_mean)
par(mfrow = c(2,2))
dem1_agg_mean = aggregate(dem1, fact = c(1,6), fun = "mean")
dem1_agg_min = aggregate(dem1, fact = 4, fun = "min")
dem1_agg_max = aggregate(dem1, fact = 4, fun = "max")
dem1_agg_sum = aggregate(dem1, fact = 4, fun = "sum")
plot(dem1_agg_mean, main = "directed aggregation")
plot(dem1_agg_min, main = "minimum")
plot(dem1_agg_max, main = "maximum")
plot(dem1_agg_sum , main = "sum")
When we want to go in the other direction, we disaggregate raster cells with the disagg() function.
Sadly, we cannot magically determine what values are the field truth for the new, smaller cells.
Instead each new cell gets the same value as it’s parent.
dem1_disagg = disagg(dem1, fact = 2)
plot(dem1_disagg)
12.1 Extracting raster values with vectors
We often encounter situations where we want to extract the values of a raster that coincide with vector data we have. For example, we might have taken soil samples along a mountain and now we want to now the altitude of the samples or we are interested in the land cover of a river catchment.
To demonstrate this, we will need some vector data sets. For the point data, we can simulate points within the bounding box of the DEMs we used before.
# - extract extend of DEM
dem.ext <- ext(dem1)
# - convert extend to bounding box usable by sf
dem.bb <- st_bbox(dem.ext)
# - create 10 random points within bounding box
random_points <- st_sample(x = st_as_sfc(dem.bb), 10)
random_points <- st_as_sf(random_points, crs = 25832)
#- interactive map
tmap_mode("view")## tmap mode set to interactive viewing
tm_shape(dem1) + tm_raster() + tm_shape(random_points) + tm_dots(size = 1)Now we have ten points that are randomly distributed across the map.
If your repeat this with you own computer your ten points will be at other positions.
Now that we have, we can extract the raster values, i.e., the altitude, at the locations of our points.
This extraction is done with extract().
The first argument to extract() is the raster raster from which we want to extract values.
The second argument is the location of points, at which we want to extract the values.
ex.points <- extract(dem1, random_points)We get a non-spatial data.frame back that holds the row number of points and the elevation at these points. We can add the altitude values back to the points like this:
random_points$altitude <- ex.points$DGM25_2905530
tm_shape(random_points) + tm_dots(col = "altitude")Now we want to look at polygons instead of points. First, we need some polygons. Here, we use river catchments from Rhineland-Palatinate.
catchments <- st_read("data/catchments.gpkg")## Reading layer `catchments' from data source `C:\Users\jonat\Documents\001_Uni\002_teaching\001_GIS\gisbook2\data\catchments.gpkg' using driver `GPKG'
## Simple feature collection with 2754 features and 55 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 4042743 ymin: 2874149 xmax: 4212789 ymax: 3094667
## Projected CRS: ETRS89-extended / LAEA Europe
Let’s have a look at the data.
tm_shape(catchments) +
tm_polygons(alpha = .5) +
tm_shape(dem_merge) +
tm_raster()